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The large-deviation method can be used to study the measurement trajectories of open quantum 
systems. For optical arrangements this formalism allows to describe the long time properties of 
the (non-equilibrium) photon counting statistics in the context of a (equilibrium) thermodynamic 
approach defined in terms of dynamical phases and transitions between them in the trajectory 
space [J. P. Garrahan and I. Lesanovsky, Phys. Rev. Lett. 104, 160601 (2010)]. In this paper, 
we study the thermodynamic approach for fluorescent systems coupled to complex reservoirs that 
induce stochastic fluctuations in their dynamical parameters. In a fast modulation limit the ther- 
modynamics corresponds to that of a Markovian two-level system. In a slow modulation limit, the 
thermodynamic properties are equivalent to those of a finite system that in an infinite-size limit 
is characterized by a first-order transition. The dynamical phases correspond to different intensity 
regimes, while the size of the system is measured by the transition rate of the bath fiuctuations. As a 
function of a dimensionless intensive variable, the first and second derivative of the thermodynamic 
potential develop an abrupt change and a narrow peak respectively. Their scaling properties are 
consistent with a double-Gaussian probability distribution of the associated extensive variable. 
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I. INTRODUCTION 

The interaction of small quantum systems with an in- 
finite set of uncontrollable degrees of freedom leads to 
non-equilibrium time irreversible evolutions. This is the 
main topic of the theory of open quantum systems [J]. 
In this context, a diverse kind of physical systems can 
be described through a quantum master equation, which 
defines the dynamics of their density matrix operators. 

Quantum optical arrangements fall in the previous cat- 
egory ■ The irreversible dynamics is induced by the 
background electromagnetic field, which leads to the nat- 
ural radiative decay of the system. The opposite mecha- 
nism, where the system becomes excited to higher energy 
levels, can be induced by the interaction with an external 
laser field. The competition between both effects pro- 
duces a continuous emission of photons. 

The detection of the radiated photons occurs at ran- 
dom times. Successive measurement realizations pro- 
vide an ensemble of trajectories, which define a stochas- 
tic (point) process [7|. Its statistics can be described 
through different approaches, such as generating opera- 
tor techniques [8j, or the quantum jump approach 
These formalisms can also be utilized in the description 
of fluorescent systems coupled to classically fluctuating 
reservoirs j9l-[llj. The interaction with the bath is mod- 
eled through a set of random processes that modify the 
characteristic parameters of the system [T2| - [l6| . In these 
cases, the main task is to relate the environment fluctu- 
ations with the photon-counting statistics [9l-[20|. 

In contrast to the previous approaches, statistical me- 
chanics provides the theoretical tools for describing sys- 
tems in thermal equilibrium [2l|. It allows to get the 
dependence of extensive thermodynamic observables as 



a function of the conjugate (intensive) variables. Phase 
transition points are defined by nonanalycities of a ther- 
modynamic potential in the parameter space. 

While both non-equilibrium (small) quantum systems 
and thermodynamic ones are described with intrinsically 
different approaches, in a recent contribution Garrahan 
and Lesanovsky related both kind of descriptions [2^ . 
The bridge between both areas is provided by the large- 
deviation (LD) theory [23j . This formalism is concerned 
with the exponential decay of probabilities correspond- 
ing to large fiuctuations in a stochastic system. It allows 
to describe its ensemble of trajectories in the same way 
as equilibrium statistical mechanics describes ensemble 
of configurations in phase space [23l - [27j . The role of the 
thermodynamic variables is played by dynamical order 
parameters or their associated conjugate fields. The exis- 
tence of "space-time" phase transitions in glassy systems 
[28j was established through this approach. 

We notice that the dynamic and physical proper- 
ties, such as subPoissonian photon statistic, photon- 
antibunching, spectrum peaks (Mollow triplet), squeezed 
noise, bistability, etc., that can be found in quantum 
optical systems are very well known Q . The break- 
through introduced in Ref. [11] is not directly related 
with those phenomena. The main idea was to apply the 
LD formalism to the measurement trajectories of simple 
open quantum systems, such as (two) three-level fluo- 
rescent systems and a micromaser. The extra physical 
aspect that can be analyzed with the LD approach is the 
asymptotic (long time) statistical properties of the mea- 
surement trajectories, such as for example the statistic 
of the number of detected photons in the case of fluores- 
cent systems or the number of atoms leaving the cavity 
in a given state for the micromaser. The LD approach. 
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by going beyond the central limit theorem [23[, allows 
to describe the asymptotic regime with a set of func- 
tions that in a statistical mechanics interpretation play 
the role of entropy and free energy. By modifying the sys- 
tem parameters and a conjugate dimensionless (dynami- 
cal) order parameter, properties such as scale invariance 
points, crossover between distinct dynamical phases, and 
an actual first-order phase transition were found in the 
thermodynamic approach j22| . As these properties are 
related with physical observables defined from the ensem- 
ble of measurement trajectories they should be detectable 
in experiments. 

The main goal of this paper is to extend the analysis 
of Ref. [ill to the case of fluorescent systems coupled to 
complex self-fluctuating reservoirs. The system-bath dy- 
namics is described through a density matrix formalism 
[iqI [20| . Our interest is to characterize the thermody- 
namic formalism associated to the photon counting prob- 
abilities as a function of the statistical properties of the 
environment fluctuations. We show that in a fast modu- 
lation limit, i.e., when the characteristic time of the bath 
fluctuations is the minor time scale of the problem, the 
thermodynamic frame reduces to that of a Markovian 
fluorescent two-level system. A scale invariance point 
is recovered for a special set of system parameter val- 
ues. On the other hand, in a slow modulation limit the 
thermodynamical properties of the measurement trajec- 
tories are equivalent to those of finite systems that in an 
infinite-size limit develop a first-order phase transition 
[29l - l32j . Here, the phases are related with different inten- 
sity regimes of the scattered radiation, while the (thermo- 
dynamic) size of the system is associated to the character- 
istic rate of the bath fluctuations. The finite-size effects 
are similar to those found in the Ising or q-state Potts 
models [3l|, [s^] ■ These properties are shown through the 
thermodynamical response functions, i.e., the first and 
second derivatives of the thermodynamic potential with 
respect to an intensive parameter. 

The paper is outlined as follows. In Sec. II, based on 
the theoretical results established in Rcf. [22], we define 
the thermodynamic approach for an arbitrary counting 
process. In Sec. Ill we define the density matrix evo- 
lution and photon counting statistics of the system of 
interest, i.e., a fluorescent (two-level) system driven by 
classical fluctuations. In Sec. IV we analyze the thermo- 
dynamic approach in the limit of fast environment fluc- 
tuations, while the case of slow fluctuations is developed 
in Sec. V. The conclusions are presented in Sec. VI. 



II. THERMODYNAMICS OF COUNTING 
PROCESSES 

Here, we define the thermodynamic formalism for 
an arbitrary counting process In the next section, 
it is build up from the photon statistics of a fluorescent 
system driven by classical fluctuations. 

A counting process is deflned by a set of trajectories. 



each one consisting in a series of consecutive events oc- 
curring at random times Q . It can be statistically char- 
acterized by a set of probabilities {-Pn(i)}^0' satisfying 
< Pn{t) < f , and the normalization 

oo 

Y.Pn{t)^l. (I) 

n=0 

Each Pn{t) is the probability of occurrence of n-events up 
to time t. From these objects, we introduce an associated 
stochastic process deflned by the probabilities 

q^{t) ^ ^P„(t)e-^", (2) 

where s is a real parameter. Consistently with the con- 
dition X]n=o'?"(^) ~ ^^"^ function Zt(s) is deflned by 

oo 

Z,(5)^^P„(t)e-^". (3) 

Hence, Zt{s) is the generating function of the original 
counting process, 

Zt(s) = ((exphsn,t(i)])){p}. (4) 

Here, nst{t) is the (stochastic) number of events up to 
time i in a given realization while {{■■■)){p} denotes 
an average over the realizations associated to the set 

From the transformation one can deduce that un- 
likely events of the counting process_[_P„(i)}^o becomes 
typical events in the s-ensemble [23|, i.e., in the set 
of realizations defined by the probabilities {qn{t)}'^=o- 
At each time t, the rare events has associated a 
thermodynamic-like structure. The consistency of this 
affirmation becomes evident after introducing the corre- 
sponding statistical objects. A thermodynamic entropy 
function St can be defined as the Shannon entropy of the 
s-ensemble, 

oo 

5* = -^g„(i)log[<z„(t)]. (5) 

ri=0 

By reading the function Zt{s) as a partition function [2l|, 
we define a "free energy function" or "grand (thermody- 
namic) potential" as 

et = -log[Ztis)]. (6) 

A "internal energy" is defined as 

oo 

{{E})^^-Y,<ln{t)log[Pn{t)], (7) 
while the "average particle number" reads 

oo 

((7V)),^5]g„(t)n. (8) 

n=0 
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Then, it is straightforward to relate the previous objects 
through the thermodynamic relation 



(9) 



In fact, in units of energy where kT — 1, k denoting the 
Boltzmann constant and T temperature, and defining a 
(dimensionless) "chemical potential" fi = — s, the previ- 
ous relation arises in the description of thermodynamical 
(equilibrium) processes carried out in open systems that 
can exchange both heat and matter with their surround- 
ings [2l|. Therefore, both the energy and particle number 
can fluctuate. Consistently with a statistical derivation 
based on maximizing entropy, {{E))^ and {{N))t can be 
read as the constraints on the average energy and par- 
ticle number respectively. Here, the average, denoted as 
((• ••))(, is defined by the set of probabilities {qn{t)}^=o- 
The thermodynamic interpretation allows us to write 
the average number (extensive variable) as the derivative 
of the (pseudo) grand potential &t with respect to the 
(pseudo) chemical potential s [2l| (intensive variable), 



(10) 



while the average of the centered quadratic fluctuations 
follows from the second derivative, 

((Aiv^)), ^ {{N^)), - mr, = -^e,. (11) 

These (two) thermodynamic relations can alternatively 
be derived by writing {{N))t and ((AiV^))t as the aver- 
age and variance of the number of events up to time t 
associated to the set of probabilities {qn{t)}^=o- 

The thermodynamic frame [Eq. ^] is parametrized by 
the time t. In a long time regime, for ergodic processes, 
it is expected that all averaged quantities (strictly, all 
cumulants) become proportional to the evaluation time 
t. Hence, the normalized asymptotic average values 



(12) 



become time independent. In this regime, the partition 
function acquires a LD form [22l| 



lini Ztis) w exp[-ie(s)]. 



(13) 



Then, the previous relations [Eqs. PH)) and ([TT])] 
maintain their validity after replacing Qt — ^ ©(s) and 
((■ ■ ■))t ^ ((• • ■ ))■ Notice that 6(s) and all normalized 
average values ((• • • )) have units of [1/time]. 

The pseudo-thermodynamic structure defined previ- 
ously is associated to the probabilities ([2]). The LD 
formalism allows to relate it with an observable prop- 
erty of the original counting process, i.e., with the long 
time behavior of the probabilities {-Pn(i) iSJLo- Consis- 
tently with Eq. (jlSp. their asymptotic structure is writ- 
ten as limt_joo ^n(i) ~ exp[— tiy9( j)]. Taking into ac- 
count that the LD function (p{N) {N = n/t) also de- 
fines the asymptotic behavior of the internal energy ([7]) 



and entropy ([5]), it can be related to the grand poten- 
tial through a Legendre-Fenchel transformation [22|, [2^ , 
(f{N) = maxs[0(s) — sN], which in turn guarantees that 
0(s) has convexity properties consistent with the ther- 
modynamical interpretation. 

From the previous relations, it becomes evident that 
the thermodynamical potential 6(s) provide an alter- 
native and complete characterization of the asymptotic 
properties of the set {Pn{t)}^=o- In consequence, possible 
phase transitions happening in the thermodynamic frame 
must to be related to strong modifications in the statis- 
tical properties of the original counting process. These 
relations and its associated theoretical frame provide an 
alternative and novel approach for analyzing measure- 
ment trajectories of single open quantum systems sub- 
jected to a continuous measurement process (2^ . 



III. FLUORESCENT SYSTEMS DRIVEN BY 
CLASSICAL FLUCTUATIONS 



In many nanoscopic optical systems, such as those 
analyzed in the context of single-molecule spectroscopy 
[9l-[ll|. the randomness of the photon emission process 
arises from both the interaction with the background 
electromagnetic field and intrinsic environment fluctua- 
tions. This last effect can be modeled by a set of noises 
that modify (modulate) the parameters of the system 
evolution [l3 - [l^ . As demonstrated in Refs. [l^llol, the 
noises modeling can be reformulated through a density 
matrix formalism. Here, we present a short derivation 
of both the underlying density matrix evolution and the 
photon counting statistics. 

The environment is defined by a set of (configurational) 
macrostates, each one leading to a different system dy- 
namics. The transitions between the bath states is de- 
scribed by a classical rate equation Q 



dPRjt) 
dt 



(14) 



where Pnit) is the probability of finding the environment 
in a given state R = 1, • • ■ Rmax, at time t. The set {(f>RR' } 
define the hopping rates. The system density matrix p{t) 
is described by a set of auxihary states {pR{t)}, each 
one representing the system dynamic for each bath state. 
Then, by writing 



(15) 



and demanding the condition 



TiipR{t)]^Pn{t), (16) 
where Tr[---] denotes a trace operation in the system 
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Hilbert space, we introduce the evolution [l9| 

- ^[Hn,PR{t)]+jnC[PB.{t)] (17) 



The last line of this equation introduce a coupling be- 
tween the auxiliary states and takes into account the en- 
vironment fluctuations. The constant define the nat- 
ural decay of the system associated to each i?-bath state. 
Consistently, the Lindblad superoperator £[•] reads [l| 



£[.] 



(7 • a' 



(18) 



where {•••}+ denotes an anticommutation operation. 
The operator cr(crt) is the lowering (raising) operator be- 
tween the system eigenstates. It is modeled through a 
two- level optical transition with eigenstates {|±)}. Then, 
CT =!-)(+! and at = |+) (_|. 

In the first line of Eq. (IT71) , Hr is the system Hamil- 
tonian associated to each bath state. It reads [19, uOI 



^{a^e~'"^* +ae+''^^*), (19) 



of the system, a "generating operator" is introduced 

e(t,s)^^°° e-V")(0, (22) 

where the extra real parameter s plays the same role as 
in Eq. ©. The partition function Eq. ^ follow as 



Zt{s)=Tv[g{t,s)] 



(23) 



From the equality p(t) = g{t,s)\s^Oi the condition 
Zt{s)\s=o = 1 is consistently satisfied. 

The conditional states p*^"^ (t) can be decomposed into 
the contributions associated to each configurational state 
of the reservoir, leading to the expression 

^(^' = Er=o PR (^) ^ E« ^«(^' (24) 

Each matrix p^'' (t) defines the state of the system under 
the condition that at time t n— photon detection events 
happened, and the environment is in the configurational 
state R. Consistently, each contribution gR{t,s) defines 
the (conditional) generating operator "given^^ that the 
reservoir is in the i?-state. Its evolution reads Il9t 



where az is the z-Pauli matrix in the basis {|±)}, and 
ujji = u}q + Slor. Therefore, Wq defines the bare transition 
frequency of the system, while Sujr are the spectral shifts 
induced by the interaction with the bath. The second 
contribution takes into account the interaction with the 
external laser field. flR is the effective Rabi frequency as- 
sociated to each reservoir state, while is the frequency 
of the external laser excitation. 

The sum structure Eq. added to the local char- 

acter of the evolution p7|) . imply the presence of strong 
non-Markovian effects in the system dynamics [l^. On 
the other hand, Eq. (IT7|) does not take into account 
light assisted processes. While these phenomena appear 
in real experimental situations [ol-fTTjl , most of the results 
developed in the next sections can be easily extended to 
include such kind of effects [19, uOl ■ 



Photon-counting statistics 

The photon detection statistics can be obtained by ex- 
pressing the system density matrix (llSp as 



n— 



(20) 



Each state p*-"-* (t) corresponds to the system state condi- 
tioned to n— photons detection events The proba- 
bility of counting rt-photons up to time t reads 



P„(i)^Tr[/")(t)]. 



(21) 



This set of probabilities can be obtained through a gener- 
ating function approach [7j . Due to the quantum nature 



dt 



'^[HR,gR{t,s)] + -iRCs[gR{t,s)] (25) 

<PrR'Gr' {t, s) - ^ <pR'RSR{t, s). 



where the superoperator Cg is given by 



^s[»] = --{crV,»}+ +e" 



(26) 



The evolution ()25|) can be solved in a Laplace domain. 
From Eq. ([M]). the partition function can always 
be written as a quotient of two polynomial functions, 
Zu{s) = flu)/g(u), where u is the Laplace variable. 
Then, by using the residues theorem, the grand potential 
Q{s) [Eq. p^ ] follows from the larger root of the equa- 
tion g{u) — 0. This property allows us to get all objects 
characterizing the s-ensemble. In the following section 
we study its associated statistic in a fast and slow mod- 
ulation limits. 

In Ref. (23 [ it was demonstrated that it is always pos- 
sible to find a density matrix p{t) whose associated mea- 
surement statistics is given by the probabilities Qnit), Eq. 
([2]). In the present context, one can affirm that the s- 
ensemble statistics can be recovered from a set of auxil- 
iary states {pR{t)}^I^Y defined as 



PRit) - 



ERT^llRPRit)]' 



1/2 



(27) 



where the set of operators {Ij^} is the left "eigenoperator" 
of the evolution ([25]) . 
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IV. FAST MODULATION LIMIT 

The thermodynamic frame associated to Eq. (|23p in 
general cannot be characterized in an analytical way. 
Nevertheless, in the limit of fast and slow environment 
fluctuations the problem becomes analytically treatable. 

For the original ensemble of trajectories [Eq. (PT|) ]. 
the fast modulation limit refers to the case in which the 
environment transitions are much faster than the pho- 
ton emission process. From Eq. (ITTI) . this condition can 
explicitly be written as 



{<I>RR'} > {Ir}- 



(28) 



The constant Iji is the intensity associated to the i?-bath 
state [diagonal contribution in (IT7|) ]. i.e., the intensity 
of a Markovian two- level system with natural decay jn, 
Rabi frequency n^j, and detuning Sr [l9l. [20j. 



IR^^ 



2f72 



(29) 



where 5r = lol — i-^r- Both, ujl and ujn are defined from 
Eq. (|19p . The inequality (|28p implies that the average 
time between two consecutive photon emissions is much 
larger than the characteristic time of the bath fluctua- 
tions. Under this condition, the fluorescent system can 
be approximated by a Markovian system whose evolution 
is deflned by the average parameters [l^, [l^ 



7 = E^«^«' ^ 

R R 



(30) 



The weights {-P^} are the stationary solution of ([M 



P]f ^ hm Pnit). 



(31) 



Therefore, the generating operator [Eq. (p4)) ] can be ap- 
proximated as 

git,s)~gMit,s), (32) 

where Quiti s) is deflned by the Markovian evolution 

^^^^^ ^~l[H,gM{t.s)]+ 7/:, [Gm (i , s)] . (33) 

Here, £s[»] follows from Eq. (|26|) while H can be read 
from Eq. (|19l) after replacing ujb, w and £7^ ^> £7. The 
sub-index M indicates the underlying Markovian approx- 
imation. 

From Eqs. (|13p and ([23| . it is possible to asso- 
ciate a thermodynamic potential Qm{s) to the operator 
^_A/(i, s). Its thermodynamics cumulants are 

m)M = ((AiV^))M = (34) 



ds 



The analytical expression for Qm{s) that can be ob- 
tained from Eq. ([55)) correspond to the larger root of a 



fourth degree polynomial. When the external laser exci- 
tation is in resonance with the (Markovian) system, i.e., 
ujl = ^, the polynomial is of third order. Thus, the ex- 
pression for Qm{s) becomes much more simple. Under 
the previous condition, from Eqs. (IT3l) and (|23)) . we get 



1 
2 



f{s) 
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Vis) 



(35) 



The auxiliary function /(s), after introducing the "fugac- 
ity" z = exp(— s) [2lj, reads 

/(s) ^ (54z7f72 + ^(54z702)2 + 27(41^2 - 72)3^ ^'^ . 

(36) 

Notice that when 7 = 217 it follows Qm = 17(1 — 
e~'^/^)^ which recovers the result presented in Ref. 
j22'|. It is characterized by the scale invariant property 
((AiV2))M/((A^))A/ = 1/3, i.e., the normalized fluctua- 
tions do not depend on s. 

In Fig. 1, we characterize the thermodynamic frame 
associated to Eqs. p7)) and (|25|) in the fast modula- 
tion limit. We assume a two-dimensional configurational 
space, i.e., the bath is characterized by only two states, 
R — A, B. We take ujfj — luq and 17/? = 17, i.e., the spec- 
tral shifts are null and the system-laser interaction is in- 
dependent of the bath states. Hence, the reservoir only 
affects the natural decay of the system, {7^}. Further- 
more, the laser excitation is assumed to be in resonance 
with the system, wl = wq. Under these conditions, we 
can approximate the grand potential as 6(s) ~ Om(s) 
[Eqs. dlSI)]. 

In Fig. 1(a) we plot the average number {{N)). In the 
inset, we show the corresponding grand thermodynamic 
potential, which is obtained as the larger root of an eighth 
order polynomial function. For the chosen parameter val- 
ues, the fast modulation limit is achieved for all values of 
s. In fact, in both cases the curves are indistinguishable 
from the analytical expressions Eqs. ([55)1 and p4p . The 
average rate [Eq. pop ] is indicated in the plots. Notice 
that for the upper (blue filled squares) curve 7/17 > 2, 
while for the lower (red filled circles) curve 7/17 < 2. 
Each value corresponds respectively to an overdamped 
[2 < 7/17 < cx)] and underdamped [0 < 7/I7 < 2] regimes 
of the Markovian dynamics, Eq. p3|) . 

In Fig. 1(b) we plot the normalized fiuctuations 
((AiV2))/((iV)). This function is almost indistinguishable 
from the fitting that follows from Eqs. ([55)1 and ([M]) . In 
the limit s +00, asymptotically all curves converge to 
1. This property is fulfilled by a Poisson process Q. In 
the limit s —00, all curves converge to 1/3, i.e., the 
value of the scale invariant regime [22l |. In the curves 
associated to the red filled circles and blue filled squares 
the parameters are the same than in Fig. 1(a). On the 
other hand, the extra curves associated to the red empty 
circles and blue empty squares correspond to a different 
set of parameter values that also are in the underdamped 
and overdamped regimes respectively. From these curves 
we deduce that, in both regimes, when 7/17 2 the tran- 
sition between 1/3 and 1 occurs for higher values of the 
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FIG. 1: (color online) (a) Average number value {{N}) as 
a function of the chemical potential s in the fast modulation 
limit. The inset show the respective grand potential Q{s). (b) 
Plot of the normalized fluctuations {{AN^)) / {{N)) . All curves 
are almost indistinguishable with the fitting defined by Eqs. 
(|34|l and ((35}. The parameters of the Hamiltonian dynamics 
are ljr = ljo and n_R = Q, while for the irreversible one read 
7A/n = 2, 7s/r2 = 3 (blue filled squares), and 7^/^ = 0.3, 
7_B/n = 0.7 (red filled circles). In (b), the parameters of 
the extra curves are ^a/^ ~ 2.5, 7s/r2 = 1.64 (blue empty 
squares), and ja/^ = 0.01, ^b/^ = 0.03 (red empty circles). 
In all cases we take (j>AB = (j>BA/^ ~ 10. The average decay 
rate, Eq. (|30p . is written in each plot. 



chemical potential, s +oo. Furthermore, we checked 
that when j/il — 2 a scale invariance property (2^ is 
recovered, {{AN'^)) / {(N)) = 1/3. The validity of this 
behavior occurs for increasing values of the rates 4>rr' . 



V. SLOW MODULATION LIMIT 

When {(j)RR'} <C {Ir}, the system is able to emits a 
huge quantity of photons before a bath transition occurs 



EE Hi. Therefore, the photon emission process de- 
velops a blinking property, i.e., the intensity randomly 
changes between the set of values {Ir}- Each change in 
the intensity regime can be associated to a bath config- 
urational transition. This limit is much more interesting 
than the previous one. Signatures of a thermodynami- 
cal phase transition can be found in this regime. First, 
we describe the emission process through a stochastic 
approximation. In a second step, we characterize the 
thermodynamic associated to the s-ensemble. 

A. Stochastic approach 

In the slow modulation limit, the realizations of the 
photon counting process [i.e., the realizations associated 
to the probabilities (HI])] can be approximated by the 
stochastic process 

nst{t)^J2 f dr5RR^,(r)'^^^. (37) 

R -^0 

As in Eq. (|4]), nst{t) is the (stochastic) number of pho- 
ton detection events up to time t. On the other hand, 
Rst{T) G (1, • • • -Rmax) IS a random process that indi- 
cates which is the state of the reservoir at time r. Then, 
the contribution defined by the discrete delta function 
^RRst{r) does not vanish only when Rst{T) = R, where it 
is equal to 1. Finally, ".||(t) is the (stochastic) number of 
photon recording events up to time t corresponding to a 
Markovian fluorescent system defined by the decay rate 
7/{, transition frequency wr, and Rabi frequency Q,r, i.e., 
the parameters associated to each reservoir state. 

Consistently with the slow modulation limit, we can 
assume that Rstir) and the set {n-|j(T)} are statistically 
independent between all them. Hence, the average num- 
ber follows from Eq. ([57)1 as 

((".tW)){P} =E / dr{SRR^^^,^)^^. (38) 

R. ^ 

With (■ • ■ ) and the overbar (• • • ) symbols we denote an 
average over the realizations of Rst{t) and n|{(i) respec- 
tively. Trivially, one can write {SRR^^(^r)) = -Pr (■''), where 
{PR{t)} are the solution of Eq. (ITil) . After assuming 

that ^(0) = and Pr{0) = Pjf [Eq. dSU], we get 

((n.t(i))){p}=E^MW- (39) 

R 

Therefore, the average number (i))){p} can be writ- 
ten as a linear combination of the averages {n||(t)}, each 
one being weighted by the stationary configurational pop- 
ulations {P"^}- Taking into account that n||(t) = iRt, 
the normalized asymptotic average value {{Anst)){p} = 
limi^oo(l/i)((Aris4(i))){p}, trivially reads 

R 
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The second cumulant, defined by 



{P}> 



(41) 



We notice tliat in addition to the linear combination 
given by the first contribution, {{An'^^{t))) ^p-^ is also pro- 
portional to the intensities {Ir}- The normalized cumu- 
lant ((An^j)) =limt^oo(lA)((An2,(i))) reads 



can be obtained in a similar way. First, from Eq. ((37l 
we write the second moment as 



{{Ani,))^p}=}_^P^Anj,- 

R 



2 ^ IrJr' Jrr' 

RR' 



{{nlmiP} = [ dT r dT'{6pR^,(r)5R,,R^,(r')) The matrix /flfl, is givcu by 

uufJo Jo 



RR 
d d 
dr dr 



^K*(TX(r')]. (42) 
The average appearing in the first line can be written as 



(WM<5i?'fl.,(r')) = P{R,t;R',t'), (43a) 
= P{R,t\R',t')Pr,{t'), (43b) 
= PiR,r\R',r')P^. (43c) 



By definition, P{R,t; R' ,t') is the joint probability for 
observing successively the bath in the states R' and R 
at times r' and r respectively (r > r'). By using the 
Markov property of the underlying bath fluctuations it 
can be expressed in terms of the conditional probability 
P{R, t\R', t') Q- As before, for simplifying the analysis, 
in the third line of the previous equation we assumed that 
the bath fluctuations begin in their stationary state. 

The second line of Eq. (|42)) define the correlation of 
the counting processes {np{T)}. Due to the statistical 
independence of these objects, when R ^ R' it follows 
n||(T)n||, (r') = ?i||(t)7i||, (t') ~ IpIpiTr'. After some 



manipulation, the second cumulant (|41l) reads 
{{AnlmiP} ^2Y,jy £dT'^^A^ir,r') (44) 
X P{R, t\R, t')P^ + 2 ^ IrIr.Jrr' (t), 

RR' 

where we have introduced the matrix of functions 

fRR'it)^ f dr r dT'[P{R,T\R',T')~P^]P^, (45) 
Jo Jo 



and the diagonal correlation A||(t, t') = np{T)np{T') — 
np{T)np {T'). The correlation time of this object is much 
smaller than the bath transition time. Hence, we can 
approximate (d/dT')An||(T, r') « S{t - t')A7i||(t, r) — 
S{t -T')An%{T), where 



Anlir) = n^^iT)n''^^{T) - n|*(r)n|*(r), 



(46) 



is the second cumulant of each process n|j (r) . From Eq. 
(HU we get the final expression 



{{Anlm^py = J2 PH^<{t) + 2 E liilR'fBii'ii)- 

(47) 



RR' 



Irr' 



lim -fRR'it), 



the intensities [Ir] follows from Eq. (|29|) and 



Anl ^ Ir 



1 - 



(7. 



2^1 



(48) 



(49) 



(50) 



where as before 6r = lvl — i^r- This expression, as well 
as that for the set {Ir}, satisfy the relations 



aew(s) 



ds 



Anl 



s=0, 



(51) 



s=0, 



where Q\,j is the grand potential associated to Eq. ([55)1 
under the replacements w — > w/f , — >■ flu, and 7 jr. 

For a two-dimensional configurational bath space, R — 
A, B, the constants can easily be written as /rr' — 
{2Srr' - l)P^P^. From Eq. (gS]), we get 



{{^nl)){P} 



R.=A,B 



PTPW{Ia-Ib? 



'JAB 



aBA) 



(52) 

where P^ = (Pab/ {<Pab + <I)ba) and P|° = (t>BA/{(f>AB + 
4>ba)- We notice that these expressions recover the results 
obtained in Ref. through a different approach. 



B. Extension to the thermodynamic approach 

The previous stochastic approach describes the s- 
ensemble in s = 0, i.e., in a slow modulation limit it can 
fit the statistical behavior dictated by the probabilities 
{Pn{t)} [Eq. dSI])] or equivalently the states {pR{t)} [Eq. 
(ITS]) ]. For s 7^ 0, it does not apply. In this case, the rel- 
evant objects are the probabilities {qn{t)} [Eq. ([2])] and 
the states {pR{t)} [Eq. (|27l)]. In the following calcula- 
tions, we assume that the slow modulation limit and the 
stochastic approach can also be applied for any value of 
the pseudo chemical potential s. The consistency of this 
ansatz relies on the results that can be obtained from it. 

By noting that ((n^t)){p} = ((iV))|s=o, from Eq. (gO]), 
we write the first s-dependent cumulant as 



(53) 



where consistently {{Np}) = ((iV))M- Here, {{N))\.j is 
defined by Eq. (p4|) after replacing w —i' up, ft ftp, 
and 7 — >■ 7_R. On the other hand, here the weights are 



P^{s)^ \imTr[pp{t)], 

t^OC 



(54) 
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where PR{t) is defined by (|27|) . 

Taking into account that ((An^j)){p} — {{AN^))\s^q, 
from Eq. ([i7|) we propose the expression 

{{AN')) ^ ^P;^(.)((A7V^^) (55) 

R 

+2j2{{NR)){{Nn,))fRR,is), 



RR' 

where ((AiV|)) = ((ATV^-^)'^ also follows from Eq. 
after replacing uj lur, fl VIr, and 7 ^ 7_r. 

The expressions ([53| and (|55l) depend on the station- 
ary populations (1541) and the generalized matrix /rr'^s). 
Its definition can read from Eq. (|49p after introducing 
in Eq. (HSl) the s-dependence of the bath populations 
[Pr — ^ Pfl(s)]. The exact expressions for these objects 
for arbitrary bath spaces [Eq. (jl4p ] are very complicated 
and do not provide an intuitive frame for understand- 
ing the thermodynamics of the s-ensemble. Therefore, 
from now on we restrict to a two-dimensional reservoir, 
R = A,B. From Eq. we write 



((AiV^)) ^ Pr{s){{ANI))+2 



R=A,B 



x{{{Na)) - {{Nb)))' 



(56) 



where (t>{s) = (1)ab{s) + 4>ba{s). 

Eqs. ([53|l and ([56|) only depends on the stationary pop- 
ulations {Pr{s)} and the generalized rate 4>{s). These 
functions can be approximated in the following way. Con- 
sistently with an equilibrium thermodynamic approach, 
the populations are expressed as 



Z{s) ■ 



-[CB+iB(s)] 



(57) 



where (r + ^r{s) is the "energy" associated to each R- 
bath state. This splitting is defined such that £,r{0) = 
0. On the other hand, the function Z{s) guarantees the 
normalization P^{s) + P^{s) = 1. Then, we can rewrite 



PTis) 

Psi^) 



{1 - tanh[£o + e(s)]}, (58a) 
{H-tanh[£o + £(s)]}, (58b) 



where eo = (Ca - Cs)/2, and e(s) = [Uis) - es(s)]/2. 
In s = 0, these expressions must to satisfy the condition 
P^(0) — P^, which imply the expression 



1 



£0 



■log 



r>oc 



poo 
^A 



1 



■log 



oba 



OAB 



(59) 



On the other hand, the dependence of e{s) is assumed to 
be linear in s, i.e., e(s) ~ sfe. The constant 5e can be 
determine from the relation between Eqs. (|53l) and ([55)) 



in s = 0, i.e., {dlds){{N))l 
some algebra, from Eq. ([SB)) . 

e{s) crL s 6e — s 



=0 = -{{^N'))l 
it follows 

Ia-Ib 



iiAB + <Pba 



After 



(60) 



Both £0 and 5e may assume positive and negative values. 

The generalized rate (j){s) [Eq. (|56p ] defines the char- 
acteristic decay time of the s-dependent bath popula- 
tions. It must to satisfies the conditions (f){s) > and 
^(0) = (0AB + 0ba)- We assume the dependence 



4){s) ~ as{s - 2sp) + {(j)AB + 4>ba), 



(61) 



where Sp is a dimensionless parameter while a > 0, as 
well as 4>{s), has units of [1/time]. As a function of s, 
the rate 0(s) reaches its minimal value at s = Sp. This 
constant is chosen as the value of s at which the second 
contribution of Eq. (I56p reaches its maximal value. It 
can be approximated by the value at which the function 
P^(s)P|°(s) = l/(2cosh[£o-h£(s)])2 is maximal. Hence, 
we take 



£0 
'6e 



'^AB 



OBA 



2{Ia - Ib) 



log 



pc 



(62) 



Finally, the coefficient a can be determine from the 
condition {d/ds){{AN''))\,=o = -((AiV3))|,^o- Alterna- 
tively, due to its large size expression, it can be considered 
as a free fitting parameter. 

The s-extension of the stochastic approach allows us to 
get closed expressions [Eqs. (|53)) and (|56)) ] for approxi- 
mating the thermodynamic approach in the slow modula- 
tion limit. In Fig. 2 we plot the average number and the 
normalized fiuctuations obtained from Eqs. (fT3)) . ([23)) 
and ([25)) . As in the previous section, the configurational 
bath space is two-dimensional, R = A,B. It only affects 
the natural decay of the system, {7^}. Hence, we take 
ujR, ~ loq, il,R = O, and ljl = wq, i.e., the spectral shifts 
are null, the system-laser interaction is independent of 
the bath states, and the laser excitation is in resonance 
with the system. 

In Fig. 2(a) we plot {{N)) as a function of s. The full 
(black) curve is the (exact) numerical solution obtained 
from Eqs. ([13]), IM]), and The grand potential Q{s) 
follows from the larger root of an eighth order polynomial 
function. The dotted (black) curves (indistinguishable) 
correspond to the fitting arising from the stochastic ap- 
proach, Eq. ([53)) . In the scale of the plot, we note that 
the average number can be approximated as 



m) 



{{Na)) 
{(Nb)) 



for s < 
for s > 



[I A >Ib)- (63) 



The superposed blue squares and red circles curves corre- 
spond to {{Na)) and {{Nb)) respectively. These contri- 
butions follow as the first derivative with respect to s of 
the grand potential 9m (s), Eq. ([35)) . under the replace- 
ment 7 — >■ 7fi. The (crude) approximation ()63p implies 
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FIG. 2: (color online) Average number value {{N}} and nor- 
malized fluctuations {{AN'^)) / {{N)) as a function of s in the 
slow modulation limit. The insets show the behaviors around 
s « 0. The full (black) curves correspond to the numerical 
solutions associated to Eqs. (|13|l . (|23|l . and (|25p . The dotted 
(black) curves are the fitting obtained from the stochastic ap- 
proach, Eqs. ((53)) and ((56|l . In (a) the blue squares and red 
circles correspond respectively to {{Na)) and {{Nb)) deter- 
mined from Eqs. p4|l and (|35|l . In (b) the blue squares and 
red circles correspond to {{ANl))/{{N)) and ((AiV|))/((iV». 
The parameters are ujn = <^o, — Q, ojl = <^o, "/a/^ ~ 2.5, 
7s/0 = 0.5, (f)AB/n = 4 X 10"^, (jisA/n = 8 x lO-"*, and 
a/n = 2.15. 



that [see Eq. 1^] P^{s) w 9{-s) and P|°(s) w 9{s), 
where 9{s) is the step function [9{s) = for s < and 
0(s) = 1 for s > 0]. As the values of the decay rates 
[7^/f2 = 2.5, ^b/^ = 0.5] correspond to the average val- 
ues of Fig. 1, the behavior of {{Na)) and {{Nb)) for any 
value of s can also be read from that plot. 

Consistently with the approximation ()63p . the behav- 
ior of {{N)) around s « seems to be discontinu- 
ous. We notice that a similar behavior arises in ther- 
modynamical first-order transitions. For example, if a 
transition is driven by temperature, the discontinuity 
in the derivative of the thermodynamic potential may 



corresponds to the difference of specific volume of two 
coexisting phases [2l|. Here, the "jump" in {{N)) is 
{{{Na)) - {{Nb)))U=o =Ia~Ib- Hence, two thermody- 
namic phases can be associated to the intensity regimes 
defined by {{Na)) and {{Nb)). 

While the average {{N)) seems to be a discontinuous 
function on a large s-scale, around the origin it is a con- 
tinuous function of s. This property is shown in the inset 
of Fig. 2(a). Even at those small scales, the stochastic 
approach, Eq. ([55)1 . provides an indistinguishable fit- 
ting (black dotted curve). In thermodynamical systems, 
finite-size effects produce a similar smoothing of the free 
energy derivative (29l - [3^ . In the present case, this rela- 
tion is established in the following section. 

In Fig. 2(b) we plot the normalized fluctuations 
{{AN^)) / {{N)) as a function of s. Consistently with the 
rough approximation (I63p . on larger s-scales we expect 
the vahdity of the approximation 



{{AN')) 



/ ((AiVD) 
l((A^I)) 



for s < 
for s > 



(I A >Ib), (64) 



where ((AA^^)) follows from Eqs. ^ and ([33]). The su- 
perposed blue squares and red circles curves correspond 
to this approximation. While they provide a very good 
fitting for \s\ > 0, around the origin the fiuctuations de- 
velops a narrow and abrupt peak (see the inset). The 
stochastic approach (black dotted line) also fits this prop- 
erty. 

The background behavior and the peak in {{AN'')) 
can be read from Eq. (j56l) . In fact, the crude approx- 
imation (I64p is indistinguishable from the contribution 

J2b=a,b PRis){{^N^))- On the other hand, the narrow 
peak is fitted by the contribution proportional to the 
product of the stationary populations P^{s)P^{s) = 
l/(2cosh[£o + e(s)])^. Therefore, the maximal value of 
the peak occurs at s = Sp, Eq. ([^ . i.e., the value of s 



the 



at which P^(s) = P^(s) = 1/2. The inset of Fig. 2(b) 
confirms this prediction. Furthermore, from Eq. ([5 
value of {{AN')) at Sp can be approximated as 



{{^N'))U 



E 

R=A,B 



Anl 



'JAB 



OB A 



y 



(65) 



[{{N))l=,^^{lA+lB)/2i 
CTr,! can be estimated as 



while the width of the peak, 



VAB 



OB A 



\Ia-Ib\ 



(66) 



Taking different values of the parameters of the evolu- 
tion ([^5)) . we have checked that in the slow modulation 
limit the position, height and width of the peak obey the 
scaling defined by Eqs. (p^ . (1551) . and (1551) respectively. 
From these expressions, one can deduce that in the limit 
i'f'AB + 'Pba) — ^ 0: the peak becomes proportional to a 
delta Dirac function. Therefore, asymptotically a first- 
order transition happens. The thermodynamic response 



10 



functions, for all values of s, are given by Eqs. (|63|) and 
i.e., the grand potential is 0(s) = Oa(s) for s < 0, 
e(s) = es(s) for s > 0, with 6(0) = 0. 



C. Finite-size effects and double-Gaussian 
approximation 

Finite-size effects in first-order transitions l29l Isoj has 
been analyzed for swtems such as the Ising [31| and q- 
state Pott models [32|. While in these systems the tran- 
sition is driven by a magnetic field or temperature, the 
thermodynamic functions have similar behaviors to those 
shown in Fig 2. The scaling of the peak in the second 
derivative of the thermodynamic potential (32j is similar 
to those of Eqs. (|62)) . ((65)) . and (IB6)) . 

From a theoretical point of view, finite-size effects 
at first-order transitions can be characterized over 
the basis of (equilibrium) Einstein fluctuation theory 
21| . which provides the probability distribution of the 
thermodynamic variable fluctuations. For example, 
for an open (thermal) system, the probability distri- 
bution V{N) of the particle number is a Gaussian 
distribution ViNj = [2TTkT{d/dn)N]-^^^ exp[-{N - 
Nf/2kT{d/dn)N], where T is the temperature, N is the 
average particle number and /i is the chemical potential. 

In the present approach, the transition is driven by the 
pseudo chemical potential s and the size of the system 
must to be inversely proportional to the rate of the bath 
fluctuations. Consistently with the Einstein fluctuation 
theory, we search for a probability distribution V{N), 
with J_^ViN)dN = 1, such that the average number 
can be obtained as 



m) = 



+ 00 



r{N)NdN, 



(67) 



while the second cumulant follows from 



((AiV')) - / ^(^)(^ - imfdN. (68) 

T(s) J-oo 

Here, the inverse of T(s) measures the "size" of the sys- 
tem [32|. In an infinite-size limit, V{N) must be a Gaus- 
sian distribution ^l|. Nevertheless, when finite-size ef- 
fects are considered in a first-order transition, one must 
to consider a superposition of Gaussian distributions, 
each one representing the coexisting phases [32| . In fact, 
the different phases are randomly explored by the system 
when its size is finite [32]. This effect is similar to the 
blinking property of the slow modulation limit [19l . [20j . 

In our problem, the coexistent phases correspond to 
the different intensity regimes defined by {{Nn)). Then, 



V{N) 



V2vr((A7Vl))T(. 



: exp 



iN-{{NA))f 



: exp 



2{{ANl))T{s\ 



(69) 



2((A7V|»T(.). 



ex. 

a 




N/Q 

FIG. 3: (color online) Probability distribution (|69|) for dif- 
ferent values of s. The vertical dotted lines correspond to 
N/Q = Ia/Q = 0.303, and N/n = Ib/Q ^ 0.222. The pa- 
rameters are the same than in Fig. 2. 



Here, {{Nr)) and ((ATV^)) follows from Eq. §^ while 
T(s) = 0(s)/2 [Eq. dSl])]. Furthermore, the weight 
of each Gaussian distribution is expressed in terms of 
the thermodynamic potential of each phase 32| . There- 
fore, the populations P^{s) are written as P^{s) ~ 
exp[—(ji — Qji{s)/(f){s)]/Z{s), where Z{s) guarantees the 



normalization P^(s) 



P^is) 



1 [compare with (|57p ] 



and 6_r(s) is the grand potential of each phase. For the 
example shown in Fig. 2, these functions can be read 
from Eq. (j35p under the replacement 7 — >■ jr. Trivially, 
the populations P^{s) can be rewritten as in Eq. ((58|) 
with Eq defined by Eq. ([59| and 



sis) 



9A(5)-9a(s) 



Ia-Ib 



+ Ois^). (70) 



This result recovers Eq. and proofs the consis- 

tency of the previous results and scaling. In fact, the 
double-Gaussian probability distribution (|69l). through 
Eqs. (p7|) and (pS)) . recovers the expressions of the ex- 
tended stochastic approach, Eqs. (|53)) and (|56| respec- 
tively. 

In Fig. 3 we plot the distribution (p^ for different 
values of s. The parameters of the underlying evolution 
are the same than in Fig. 2. Near of the transition, s « 
Sp w 0, the probability distribution is a double-Gaussian 
one. Consistently, the higher peaks are centered around 
the values iV w /a (s < 0) and N ^ Ib (s > 0). For 
\s\ ^ Sp, the distribution has only one peak, which is 
centered around s w {{Nr}). 

The double-Gaussian approximation, and consistently 
the stochastic approach, can be extended beyond the slow 
modulation limit. Nevertheless, parameters such as the 
position Sp and width (jp of the peak must be taken as 
free parameters. For example, a reasonable fitting is ob- 
tained after replacing the polynomial function (j6ip with 
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an hyperbolic one and introducing a non-linear function 
e(s), both of them defined with extra free parameters. 

VI. SUMMARY AND CONCLUSIONS 

The (non-equilibrium) ensemble of measurement real- 
izations of an open quantum system can be analyzed with 
the LD formalism. For fluorescent systems under a direct 
photon detection scheme, this approach allow to describe 
the asymptotic behavior of the photon counting proba- 
bilities through a thermodynamic-like formalism [22l | . In 
this paper we have studied the thermodynamic approach 
associated to a fluorescent system coupled to a complex 
self-fluctuating environment able to modify the charac- 
teristic parameters of the system evolution. 

The statistical mechanics underlying the thermody- 
namic frame is defined by a set of auxiliary probabili- 
ties whose characteristic events are the unlikely ones of 
the photon counting realizations, Eq. ([2]). A free energy 
function, Eq. through a Legendre transformation, 

defines the long time behavior of the photon counting 
probabilities. Here, its functional form follows from the 
trace of a generating function operator, Eq. whose 
evolution, Eq. (l25t . takes into account the parameter 
fluctuations induced by the environment. 

In a fast modulation limit, i.e., when the characteristic 
time of the environment fluctuations is much smaller than 
the time between consecutive photon emissions, the ther- 
modynamic frame can be well approximated with that 
corresponding to a Markovian fluorescent system. Its 
evolution is defined by a set of parameters that follows 
from an average weighted by the stationary populations 
of each bath state, Eq. (15(11) . When the bath only affects 
the natural decay of the system and the external laser 
excitation is in resonance with the system, the thermo- 
dynamic potential can be approximated by a simple an- 
alytical expression, Eq. psp . In this limit, the response 
functions do not display any property related to a phase 
transition. Nevertheless, the normalized fluctuations al- 
ways interpolate between a scale invariant regime and a 
Poissonian one [Fig. 1]. 

In the slow modulation limit, the fluorescent signal is 
characterized by a blinking phenomenon. The scattered 
intensity randomly changes between a set of values asso- 
ciated to each bath state, Eq. (I^H]) . The photon counting 
process can be approximated by the product of two kind 
of statistically independent stochastic variables, one re- 
lated to the quantum photon emission process and the 
other to the bath fluctuations, Eq. ([57)) . In the ther- 



modynamic frame, each intensity regime can be read as 
a different thermodynamic phase. A natural extension 
of the stochastic approach provides the basis for charac- 
terizing its statistical properties. After imposing some 
consistency relations, the average number is written as 
a linear combination of the values corresponding to each 
phase, Eq. ([55]) . The fluctuations around the average 
number can be approximated in a similar way, Eq. (j56p . 

The behavior of the average number and the centered 
fluctuations is similar to that found in flnite-size systems 
near a first-order phase transition. This is the main result 
of this contribution. Instead of a discontinuity in the first 
derivative of the thermodynamic potential, an abrupt but 
continuous change in its slope is observed. Furthermore, 
the second derivative, instead of a delta Dirac contribu- 
tion, displays a narrow peak [Fig. 2]. These effects are 
controlled by the size of the system, which is proportional 
to the transition rate between the bath states. The lo- 
cation, height and width of the peak obey the scaling 
properties obtained from the stochastic approach, i.e., 
Eqs. ([nil), (|55|) . and ([M)) respectively. The finite-size ef- 
fects can also be obtained from a generalization of the 
Einstein's fluctuations theory. The probability distribu- 
tion of the fluctuations follows from a double-Gaussian 
distribution [Eq. (p^ ]. each contribution being related 
to each coexisting phase [Fig. 3]. 

From our results, we conclude that whenever a (pho- 
ton) counting process has an underlying blinking prop- 
erty, the thermodynamic approach is characterized by 
flnite-size effects corresponding to a flrst-order transi- 
tion. Therefore, the studied phenomena, for example, 
must also appear when the blinking properties depend 
on the external laser excitation, i.e., for light assisted 
processes [l^ [l^j . These and previous results Q con- 
flrm that diverse thermodynamical properties of many 
body (equilibrium) systems are also present in the statis- 
tical properties of (non-equilibrium) quantum measure- 
ment trajectories. This mapping raise up fundamental 
physical questions such as the possibility of simulating 
complex dynamics with simple open quantum systems 
subjected to a continuous measurement process. 
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